Technical  Research  Report 


Coordinated  Orbit  Transfer  for  Satellite  Clusters 


by  Fumin  Zhang,  P.S.  Krishnaprasad 


TR  2002-35 


I  SR  develops,  applies  and  teaches  advanced  methodologies  of  design  and  analysis  to  solve  complex,  hierarchical, 
heterogeneous  and  dynamic  problems  of  engineering  technology  and  systems  for  industry  and  government. 

ISR  is  a  permanent  institute  of  the  University  of  Maryland,  within  the  Glenn  L.  Martin  Institute  of  Technol¬ 
ogy/A.  James  Clark  School  of  Engineering.  It  is  a  National  Science  Foundation  Engineering  Research  Center. 

Web  site  http://www.isr.umd.edu 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

2QQ2  2.  REPORT  TYPE 

3.  DATES  COVERED 

4.  TITLE  AND  SUBTITLE 

Coordinated  Orbit  Transfer  for  Satellite  Clusters 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Office  of  Scientific  Research, 875  North  Randolph 

Street, Arlington, VA, 22203-1768 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

see  report 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF:  17.  LIMITATION  OF 

18.  NUMBER  19a.  NAME  OF 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE 

unclassified  unclassified  unclassified 

20 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


Coordinated  Orbit  Transfer  for  Satellite  Clusters* 


Fumin  Zhang  and  P.S.  Krishnaprasad 

Institute  for  Systems  Research  & 
Department  of  Electrical  and  Computer  Engineering 
University  of  Maryland  at  College  Park 
College  Park,  MD  20742 

January  2,  2003 


Abstract 

We  propose  a  control  law  which  allows  a  satellite  formation  to 
achieve  orbit  transfer.  During  the  transfer,  the  formation  can  be  either 
maintained  or  modified  to  a  desired  one.  Based  on  the  orbit  transfer 
control  law  proposed  by  Chang,  Chichka  and  Marsden  for  single  satel¬ 
lite,  we  add  coupling  terms  to  the  summation  of  Lyapunov  functions 
for  single  satellites.  These  terms  are  functions  of  the  difference  between 
the  mean  anomalies  (or  perigee  passing  times)  of  formation  members. 
The  asymptotic  stability  of  the  desired  formation  in  desired  orbits  is 
proved. 


1  Introduction 

This  paper  is  concerned  with  the  problem  of  achieving  a  satellite  formation 
near  a  designated  elliptic  orbit.  For  orbits  near  the  earth,  one  can  use  a 
space  shuttle  to  place  satellites  into  specified  relative  positions.  What  we 

‘This  research  was  supported  in  part  by  the  National  Aeronautics  and  Space  Adminis¬ 
tration  under  NASA-GSFC  Grant  No.  NAG5-10819,  by  the  Air  Force  Office  of  Scientific 
Research  under  AFOSR  Grant  No.  F49620-01-0415,  by  the  Army  Research  Office  under 
ODDR&E  MURI97  Program  Grant  No.  DAAG55-97- 1-0114  to  the  Center  for  Dynam¬ 
ics  and  Control  of  Smart  Structures  (through  Harvard  University),  and  under  ODDR&E 
MUR101  Program  Grant  No.  DAAD19-01-1-0465  to  the  Center  for  Communicating  Net¬ 
worked  Control  Systems  (through  Boston  University). 
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want  to  consider  here  is  the  case  when  members  of  the  formation  have  been 
placed  relatively  far  apart.  They  have  to  use  their  on-board  thrusters  to  get 
to  the  desired  orbits  to  form  the  desired  formation.  A  similar  case  is  when 
the  whole  formation  has  to  be  restructured  for  mission-related  reasons.  We 
want  the  formation  to  be  maintained  to  some  extent  during  the  transfer  and 
be  re-established  after  the  transfer. 

Our  approach  is  to  use  Lyapunov  functions  to  design  the  control  laws 
for  orbit  transfer.  The  Lyapunov  function  will  achieve  a  local  minimum 
when  correct  orbit  and  formation  are  reached.  In  [4],  a  Lyapunov  function 
is  expressed  as  a  quadratic  function  of  the  differences  of  orbital  elements 
between  current  orbit  and  the  destination  orbit.  However,  the  convergence 
of  the  associated  control  law  is  not  proved. 

In  this  paper,  we  develop  a  new  control  algorithm  and  give  a  proof  of 
convergence.  This  algorithm  is  based  on  a  Lyapunov  function  on  the  shape 
space  of  elliptic  orbits  in  our  previous  work  [5],  the  work  of  Chang,  Chichka 
and  Marsden  [1]  and  a  result  of  Cushman  and  Bates  [2] .  However,  the  most 
significant  extension  is  the  addition  of  a  coupling  term  which  is  a  function 
of  the  difference  between  the  osculating  perigee  passing  times. 

In  section  2,  we  develop  formulas  used  in  the  proofs  of  our  theorems.  In 
section  3,  a  brief  summary  of  results  in  [1]  and  [5]  are  given. We  introduce  the 
definition  of  periodic  satellite  formations  in  section  4.  Our  main  results  and 
proofs  about  orbit  transfer  of  periodic  formations  are  presented  in  section 
5.  Simulation  results  are  shown  in  section  6. 

2  Preparations 

If  the  mass  of  a  satellite  is  small  compared  to  the  mass  of  the  earth,  the 
Kepler  two  body  problem  can  be  approximated  by  a  one  center  problem  as: 

mq  =  — V  Vq  +  u  (1) 

where  q  £  7L:!  is  the  position  vector  of  the  satellite  relative  to  the  center  of 
the  earth,  m  is  the  mass  of  the  satellite,  Vq  is  the  gravitational  potential  of 
the  earth,  u  is  the  control  force  plus  other  disturbances.  Without  considering 
higher  order  terms,  Vq  takes  the  form 

VG  =  (2) 

II  Q  II 

Let  p  =  mq  be  the  momentum  vector  of  the  satellite.  For  simplicity  we 
assume  that  all  the  satellites  considered  in  this  paper  have  unit 
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mass. 

Let  us  make  the  following  definitions: 


m 

m 

e(t) 

a(t) 


cos(E(t)) 

M(t) 


q{t)  X  p(t) 
p(t )  X  l(t)  -  p 


g(*) 

?(*) 


A(t) 


h(t)2 


Ml-e(t)2) 
e(f)V  a(f); 

-E(t)  —  e(t)sin(E(t )) 


(3) 


where  /i(t)  =  ||  l(t)  ||  and  r(t)  =  ||  q(t)  ||.  These  formulas  can  be  found  in 
textbooks  on  celestial  mechanics  [3].  I  is  the  angular  momentum  vector.  A 
is  called  the  Laplace  vector.  They  are  conserved  if  u(t)  =  0.  a  is  the  length 
of  the  semi-major  axis  and  e  is  the  eccentricity.  E  is  the  eccentric  anomaly 
and  M  is  the  mean  anomaly.  The  last  equation  is  Kepler’s  equation.  When 
e(t)  =  0,  the  eccentric  anomaly  E(t)  is  defined  to  be  M(t).  For  now,  we  will 
assume  that  e(t)  /  0. 

When  u(t)  /  0,  the  quantities  defined  above  are  called  osculating  ele¬ 
ments.  Notice  that  these  formulas  are  valid  for  all  t  and  all  the  elements  are 
differentiable  on  TZ3  x  TZ3  —  {0}.  So  we  can  take  derivative  on  both  sides  of 
equations  (3).  Take  l(t)  for  example,  we  have 


m 


dl  dl 

d^q+YPP 

dl  .  d l  .  

d-qq+dp^VVG  +  u) 
dl 

l\u=  0  +  g^u{t) 

dl  ,  . 

a~Puit) 


(4) 


Notice  that  we  used  the  property  that  l  is  conserved  when  u(t)  =  0.  Simi¬ 
larly,  we  have 
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Next, 


e{t)  = 


A- A 

OTUI 

1  cM, 

H  dp 


A  dA 

HOT  '  ~dvu[t) 


=  -(“ST  )TA-u(t) 


where  A  is  the  unit  vector  along  the  direction  of  A.  We  also  have 

:  z  \  l  *  l  1  /  01 ,7^  ,  , 

( ’ = pi = sW  ■ n(t) 


(6) 

(7) 


d(t)  = 


HI  -  e2) 


(2hh  +  2  paee) 


2  HCt, 

:(o-)  l-u(t) 


p(l  —  e2)  Hp 

2ae  ,dA,r  -  .  . 

+W^F)(^)  A'u(t) 

Since  r2  =  q  ■  q ,  we  have  rf  =  q  ■  p.  Hence 

h2  =  l  ■  l  =  (q  x  p)  ■  (q  x  p) 

=  {q  ■  q){p  ■  p)  -  (q  ■  p)(q  ■  p) 

=  r2p( - )  —  r2r2 

r  a 


On  the  other  hand 

Thus 


h2  =  pa(  1  —  e2) 

■  2  _  pa(l-e2)  2  1 


+  M - ) 

r  a 


Apply  r  =  a(l  —  ecos(-E)),  after  simplification  we  have 

_ e  sin(E) 

r  =  y/Jfd - 

r 

Taking  derivative  on  both  sides  of  r  =  o(l  —  ecos(E )),  we  have 
f  =  a(l  —  e  cos(E))  +  a(e  cos(E)  —  e  sin(E)E) 

Thus 

r 


E  = 


r  acos(E) 

-a  H - 


aesin(E)  a2esin(E)  aesin(E) 


(8) 


(9) 

(10) 

(11) 


(12) 


(13) 
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cos(E)  .  r 

esin(E)  a2esin(E) 


(14) 


On  the  other  hand,  taking  derivative  on  both  sides  of  M  =  E  —  esin(E), 
we  get 


M  =  E  —  esin(E)  —  ecos(E)E 

=  (1  —  ecos{E))E  —  esin(E)  (15) 


Combining  this  equation  with  equation  (14),(8)and  (6),  we  have 


where  n 


M  = 


cos(E)(  1  —  ecos(E)) 


e  sin(E) 


(1  -  ecos(E))J^~  +  ( 

V  a  r 

■  f  tpw  ’  r(1  -ecos(E))  . 

v  ”  a2e  sin(E) 

JT  cos(E)  —  e.  (1  —  ecos(E))2  . 

H - •  7  e - •  ,  T-'\ — a 

ae  sin(E) 


aJ  esin(E) 


=  n  +  v{^)TA-u{t)  +  £{^)Tl-u(t) 


2tt/T  and 

v(l,A,E ) 


2(1  —  ecos(E))2 
/j,ae(  1  —  e2)sin(E) 
cos(E)  —  e  2(1  —  ecos(E))2 

Hesin{E )  fi(\  —  e2)sin(E) 


(16) 


(17) 


Notice  that  £  and  i)  will  be  oo  if  sin(E)  =  0.  Physically  it  means  that  when 
the  satellite  is  passing  perigee  or  apogee,  the  control  can  cause  a  sudden 
jump  of  the  mean  and  eccentric  anomalies.  In  order  to  prevent  this  from 
happening  in  our  control  laws,  we  will  turn  off  the  control  when  sin(E)  =  0. 


3  Shape  space  and  orbit  transfer  of  single  satellite 

For  a  single  satellite  on  an  elliptic  orbit,  the  set  D  of  ordered  pairs  (l,  A)  is 
a  subset  of  7 £3  x  7 £3  with  Euclidean  norm, 

D  =  {(/,  A)  €  R3  x  R3  |  A  ■  l  =  0,  l  /  0,  ||  A  ||  <  m2n}  (18) 

Let  W  =  Tj'm  ||  ||  +  Vq  be  the  total  energy  of  the  satellite.  Let  P  be  the 
set  of  all  pairs  ( q,p )  with  Euclidean  norm.  Then  we  can  define  a  set  £e  as 

Se  =  {(q,p)  G  P\  W{q,p)  <  0,1  /  0}  (19) 

By  the  definition  of  l  and  A,  we  have  a  mapping  ir  :  T,e  — >  D,  ( q,p )  (/,  A). 
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Theorem  3.1  (Chang-Chichka-Marsden)  [1]  The  following  hold: 

1.  Ee  is  the  union  of  all  elliptic  Keplerian  orbits. 

2.  7r(Xe)  =  D  and  £e  =  tt~1(D). 

3.  The  fiber  it-1  ((l,  A))  is  a  unique  (oriented)  elliptic  Keplerian  orbit  for 
each  (l,  A)  £  D.  (see  also,  [2],  page  58) 

The  mapping  n  is  a  continuous  mapping  because  (l,  A)  are  continuous  with 
respect  to  (q,p). 

Corollary  3.2  7r_1(lv)  is  compact  for  any  compact  set  I\  C  D  (c.f.  [5]) 

To  control  the  orbit  transfer  of  a  single  satellite,  one  considers  a  Lyapunov 
function  from  [1] 

V(q,p)  =  ^(\\l-ld\\2+\\A-Ad\\2)  (20) 

where  (ld.  Ad)  is  the  pair  of  the  angular  momentum  vector  and  Laplace 
vector  of  the  target  elliptic(circular)  orbit.  The  derivative  of  V  along  the 
integral  curve  of  the  system  is 

V  =  {l-ld)4  +  {A-Ad)-A 

n  1  \  dZ  ,  dA 

=  Q-ld)-QjU  +  (A-Ad)-—u 

=  [(^  — 

lx(A  —  Ad)  +  ((A  -  Ad)xp)xq\  ■  u  (21) 

If  we  let  the  control  to  be 

u  =  —  [(Z  —  ld)  xq  +  lx  (A  —  Ad)  +  ((^4  —  Ad)  xp)xq]  (22) 

then  V  <  0  along  the  trajectory  of  the  closed  loop  system.  The  following 
lemma  is  proved  in  [5] . 


Lemma  3.3  Suppose  a  single  satellite  has  external  control  u(t )  =  0.  Let  x, 
y  be  time  invariant  unknown  vectors.  Suppose  ( q{t),p(t ))  £  Ee,  the  solution 
of  equation 

xxq  +  Ixy  +  (yxp)xq  =  0  (23) 


is 


x  =  aA  y  =  al 


(24) 


For  some  a  £  R 
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We  will  give  a  proof  of  the  following  theorem  by  Chang,  Chichka  and 
Marsden  [1]: 

Theorem  3.4  (Chang-Chichka-Marsden)  There  exists  c  >  0  such  that  if 
V{qo,Po)  <  c,  by  applying  the  control  law  as  in  equation  (22),  the  trajectory 
of  the  closed  loop  system  starting  at  (qo,Po)  will  asymptotically  converge  to 
the  target  orbit  7r_1((/^,  Ad)) 

Proof  Let 

Q  =  {(l  A)  g  1Z3  x  TZ3\l  -A  =  0}  (25) 

Recall  that  the  set  of  (l,  A)  pairs  for  elliptic  orbits  is 

D  =  {( l ,  A)  G  1Z3  x  fZ3\l  •  A  =  0,  l  ±  0,  ||  A  ||  <  m2ju}  (26) 

It  is  easy  to  see  that  Q  is  a  closed  subset  of  TZ3  x  7Z3  and  D  is  a  subset  of 
Q  with  nonempty  interior.  The  set 

B  =  {(/,  A)  €  n3  x  n3\V  <  c}  (27) 

is  a  closed  ball  in  TZ 3  x  7L3.This  tells  us  that  the  set  B  n  Q  is  a  compact 
subset  of  TZ3  x  TZ3.  If  c  small  enough,  we  have  B  n  Q  C  H,then  by  corollary 
3.2,  the  set  tt~1(B  n  Q)  will  be  a  compact  subset  of  Xe. 

The  condition  V(qo,po)  <  c  tells  us  that  (qo,Po)  £  tt^1(B  n  Q).  By 
applying  LaSalle’s  invariance  principle, we  conclude  that  the  trajectory  of 
the  closed  loop  system  will  converge  to  the  maximal  invariant  set  within  the 
subset  of  B  n  Q  where  V  =  0.  In  [5],  we  showed  that  the  maximal  invariant 
set  is  the  set  where  u(t)  =  0  for  all  t. 

According  to  lemma  3.3,  the  solution  to 

(l  —  Id)  xq  +  lx  {A  —  Ad)  +  ((A  —  Ad)  xpxq)  =  0  (28) 


is 


l  —  Id  =  olA 


A  —  Ad  =  al 

Since  Id  ■  Ad  =  0  and  l  ■  A  =  0,  we  have 

(29) 

a(im|2  +  M||2)  =  0 

(30) 

If  we  further  assumed  that 
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The  closed  loop  system  will  not  reach  the  set  where  l  =  0  and  A  =  0.  Thus 
the  only  possibility  for  equation  (30)  to  be  true  is  to  have  a  =  0.  So  the 
solution  is 


l  =  Id 

A  =  Ad  (32) 


The  authors  of  [1]  do  not  give  an  explicit  upper  bound  for  the  value  of  c. 
This  upper  bound  can  be  calculated  by  solving  a  constrained  maximization 
problem  as  below: 

sup{^( ||  l  -  ld  ||2  +  ||  A  -  Ad  ||2)}  (33) 

under  constraints 

(l,  A)  e  Q,  1/0  and  ||  A  ||  <  m2/i  (34) 

First,  we  need  to  calculate  the  supremum  of  the  unconstrained  maxi¬ 
mization  problem.  It  is  easy  to  see  that  this  value  is  oo.  Then,  we  need 
to  calculate  the  minimum  value  subject  to  l  =  0.  The  result  is  0.5  ||  ld  || 2 
achieved  when  A  =  Ad.  We  should  also  calculate  the  minimum  value  subject 
to  ||  A  ||  =  ?n2/i,  by  applying  the  Lagrange  multiplier  method  we  found  this 
value  to  be  0.5 (m2fr  —  ||  Ad  ||)2.  Hence,  for  theorem  3.4  to  be  valid,  we  must 
have 

c  <  min{-  ||  ld  ||2  ,  -{m2^  -  ||  Ad  ||)2}  (35) 

which  guarantees  (31). 


4  Periodic  formation 

Suppose  we  have  a  formation  consisting  m  satellites.  Let  Oj  denote  the  orbit 
of  the  jth  satellite.  The  orbit  Oj  and  the  position  of  the  satellite  on  Oj  can 
be  described  either  by  the  six  orbital  elements  ( a  j ,  ej ,  I j  Aij,  ujj ,  Tj )  or  by 
( lj,Aj,Tj )  where  a  denotes  the  length  of  semi-major  axis,  e  is  the  eccentricity, 
I  is  the  inclination,  H  is  the  longitude  of  the  ascending  node,  cu  is  the 
argument  of  the  perigee  and  r  is  the  perigee  passing  time.  If  Oj  is  elliptic  or 
circular,  except  for  some  singular  cases  when  some  of  the  orbital  elements  are 
not  well  defined,  theorem  3.1  tells  us  that  the  two  expressions  are  equivalent. 
Among  all  the  possible  formations,  we  are  interested  in  formations  with 
periodic  shape  changes.  Hence  we  have  the  following  definition. 
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Definition  4.1  A  formation  is  periodic  when  the  orbits  Oj  satisfies  aj  = 
a  >  0  for  all  j  =  1,  2,  ...m 

This  definition  is  valid  since  all  the  satellites  in  a  periodic  formation  will 
have  the  same  orbital  period 

T  _  2tt \fafi 

Thus  although  the  shape  of  the  formation  is  varying,  it  is  varying  periodi¬ 
cally. 

On  the  other  hand,  given  a  set  of  orbits  with  the  identical  length  of 
semi-major  axis,  there  are  infinitely  many  possible  periodic  formations.  We 
need  to  specify  their  relative  positions  on  the  orbits  to  determine  a  specific 
formation.  The  difficulty  is  that  the  relative  positions  are  complicated  func¬ 
tions  of  time.  However,  the  differences  (Tj  —  Tj)  are  constants.  If  we  define 
the  mean  anomaly  as 

M.i  =  nfit  -  Ti )  (37) 

where  n*  =  27 r/T  ,  ( Mt  —  Mj)  are  constants.  By  specifying  the  values  of 
(r,;  —  Tj)  or  (Mi  —  Mj)  for  all  i  and  j,  a  periodic  formation  can  be  uniquely 
determined. 


5  Control  laws  for  orbit  transfer  of  satellite  for¬ 
mations 

To  set  up  a  periodic  formation  of  two  satellites,  one  can  control  each  satellite 
separately  to  transfer  to  its  target  orbit.  However,  this  will  not  assure 
the  correct  values  of  (t;  —  Tj)  or  (Mi  —  Mj).  In  order  to  do  that,  extra 
terms  involving  (r,  —  Tj)  or  (Mj  —  Mj)  should  be  added  in  the  summation 
of  the  Lyapunov  functions  for  single  satellites.  This  extension  will  result  in 
a  cooperative  orbit  transfer  of  multiple  satellites. 

We  introduce  a  variable  Tj  which  is  defined  as 
Tj  =  J  — |- Mj ,  if  Ei  €  (— e,  7 r  +  e) 

V  ad 

Tj  =  J^(2tt-  Mi),if  Ei  €  (7r-e,27T  +  e)  (38) 

V  ad 

for  i  =  1,2,  where  a d  is  the  common  length  of  the  semi-major  axes  for  the 
destination  orbits. 
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Here,  the  trouble  of  using  different  expressions  for  the  cases  E)  G  (— e,  7T+ 
e)  and  Ej  G  (tt  —  e,  2ir  +  e)  is  caused  by  the  fact  that  Ei  £  S1  a  circle.  Two 
coordinate  charts  are  required  on  Sl.  Here  we  pick  the  charts  to  be 


'ipi  :  (-e,  7T  +  e)  — »  (-e,  7 r  +  e)  s.t.  Ei  Ei 
V>2  :  (tt  -  e,  27t  +  e)  — >  (— e,  7r  +  e)  s.t.Ei  (27r  -  E j)  (39) 


Here,  the  value  of  e  is  chosen  so  that  the  two  satellites  will  always  be  in 
the  same  chart.  Because  in  a  satellite  formation  the  angular  separations 
between  satellites  are  usually  small,  the  value  of  e  is  small. 

For  Ei  G  (— e,  7r  +  e),  we  have 


r,  =  J^Mi  +  1  J^OiMi 

2\ad 


dpi 


dh 


+C(^>  A-i,  Ei){—  )T U  •  Ui(t ) 
opt 


(40) 


where 


c  = 


l%Z(li,Ai,Ei)+3  la‘ 
a° 


2V«dMl-ef) 


Mi 


I  (1  A  7— I  \  |  ^  2 CLiCi 

p  -  + 


(41) 


For  Ei  G  (n  —  e,  27r  +  e), 


T*  =  -J^Mi  +  lMa^n-Mi) 

ad  z\l  ad 

M  +p(^,A,^)(^)tA-^W 

Opi 

f)1 

+C(luAi,Ei)(-^)Tli-Ui(t) 
opi 


(42) 


where 
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i  a :  .  .  3 

P  =  ~\j  73  V(h,  A'ii  Ei)  + 


2  a,  a 


{2tt  -  Mi) 


(43) 


Notice  that  we  have  terms  that  explicitly  contain  M, .  If  not  handled 
well, these  terms  will  cause  discontinuity  in  our  control  algorithm  when  the 
satellites  enter  a  new  chart.  The  reason  for  us  to  pick  the  particular  charts 
(V’i,Y’2)  is  to  reduce  the  discontinuities  in  the  derivatives  of  T*  caused  by 
changing  charts. 

We  will  design  a  Lyapunov  function  on  the  phase  space  of  the  two  satel¬ 
lites.  This  one  function  will  have  different  expressions  in  different  charts. 
The  Lyapunov  function  is 

V  =  V1  +  V2+4sin(Tl~J2~</))2 
if  Ei  €  (-e,  7T  +  e) 

v  =  Fi  +  y2  +  4  «(Tl"^2  +  ^)2 

if  Ei  G  (7T  -  e,  27 r  +  e)  (44) 

where 

Vi  =  2  (II  ^  —  ^dl  II2  +  II  —  Adi  ||2) 

V2  =  -(||  l2  ~  ld2  ||2  +  ||  A2  —  Ad2  ||2)  (45) 

Here,  and  (J-ds?  Ad2)  specify  the  orbits  in  a  two-satellite  periodic 

formation  and  specifies  the  desired  (Mi  —  M2)  on  these  orbits. 

We  can  calculate  the  derivative  of  V  as 

V  =  V\  +  V2  +  sin(Tl  ~  ^  T  </>)(f  x  -  f  2)  (46) 

The  choice  of  —  or  +  depends  on  the  value  of  E,  as  in  the  definition  of  V. 

By  the  calculations  performed  in  the  single  satellite  case, 

k*  —  \(Ji  ldi)^Qi  T  ^’X(j4j  A-di) 

+((Ai  -  Adi)xpi)xqi ]  ■  m  (47) 

for  i  =  1,2.  Thus 

V  =  [(h  ~  Idi  +  Ci^(Tl~I2=F<^)^i)xgi+ 

h  x (Ar  -  Adi.  +  p1sin( Tl~ J2^)li)+ 
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(04i  -  Adl  +  pisin(Tl  z2T(t,)A1)xp1)xq1\  ■  m 
+[(h  ~  ld2  -  (2sin(r  1~22T<?!>)/2)xff2+ 
hx(A2  —  Ad2  —  p2sin(— — J2=F<^)^42)+ 

((A2  -  Ad 2  -  P2sin{Tl~12^<t>)A2)xp2)xq2]  ■  u2  (48) 

In  order  to  get  V  <  0,  we  let 

u\  =  —  sin2 {E\)[{li  -  ld i  +  (isin(Tl~22T<l,)li)xqi+ 

h x (Ai  -  Ad i  +  p1sin(Tl~22:f<l,)Ai)+ 

((Ai  -  Ad i  +  pisin(Tl~12T<t>)Ai)  xpi)xqi] 
u2  =  -sin2(E2)[(l2  ~  ld.2  ~  (2sin(r^2^)l2)xq2+ 
l2x(A2  -  Ad2  -  p2sin(Tl~^2^(t,)A2)+ 

((A2  -  Ad2  -  p2sin(Tl~12T<t>)A2)  xp2)xq2}  (49) 

Notice  that  the  factors  sin2(Ei)  cancel  the  term  sin(Ef)  in  the  denom¬ 
inators  of  Q  and  pi.  This  will  result  in  a  continuous  control  law  which  will 
be  0  when  £)  =  0, 7 r. 

Let  z  =  {qi,Pi,q2,P2)-  We  now  proceed  to  find  the  initial  condition 
zo  =  (<7i(0),pi(0), <72(2), p2(0))  for  2  s.t.  the  set 

SM  =  {z\V(z)  <  V(z0)}  (50) 

is  a  compact  subset  of  £ei  x  Se2  —  {z\A\  =  0  or  A2  =  0}.  This  is  a  necessary 
step  because  we  want  to  apply  LaSalle’s  invariance  principle  to  prove  our 
main  result. 


Lemma  5.1  Let 

where 

.  A 

Ci  =  mm\  - 
l2 

for  i  =  1,2.  Then  the  set 


c  <  min{ci,  C2} 


Adi  ii24ii^ii24(m 


Adi  II)2} 


(51) 

(52) 


SM  =  {z\V(z)  <c]  (53) 

is  a  compact  subset  of  Eei  x  Se2  —  {z\A\  =  0  or  A2  =  0}.  Proof  The  first 
observation  is  that  the  set 


<51  =  {(Qi,Pi)\V(qi,pi)  <  c*,c*  <  ci} 


(54) 
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is  a  subset  of  £ei  i.e.  Si  n  £ei  =  Si. 

In  fact,  ci  is  the  supremum  of  V\ {q\ , p\ )  on  the  set  £ei  —  { (<?i , Pi ) | ^4i  =  0}. 
To  see  this,  we  solve  a  constrained  maximization  problem  as  below: 

sup{^( ||  h  -  ldi  ||2  +\\Ai-  Adi  ||2)}  (55) 

under  the  constraints 

A\  ■  h  =  0  l\  0  A\  /  0  and  ||  A\  ||  <  p  (56) 

First,  we  need  to  calculate  the  supremum  of  the  unconstrained  maximiza¬ 
tion  problem.  It  is  easy  to  see  that  this  value  is  oo.  Then,  we  need  to  cal¬ 
culate  the  minimum  value  subject  to  l\  =  0.  The  result  is  \  ||  Id  ||2  achieved 
when  A  =  Ad-  Similarly,  the  minimum  value  subject  to  A\  =  0  is  \  ||  Ad  ||2. 
We  should  also  calculate  the  minimum  value  subject  to  ||  A  ||  =  //.  By  apply¬ 
ing  the  Lagrange  multiplier  method  we  found  this  value  to  be  —  ||  Ad  ||)2. 
Thus 

ci  =  min{-  ||  Adl  ||2  ,  -  ||  ldl  ||2  ,  -{p  -  ||  Adl  ||)2}  (57) 

is  the  supremum  of  V\  on  the  set  £ei  —  {(gi,pi)|^4i  =  0}.  We  proved  that 

51  C  £ei  -  {(qi,pi)\Ai  =  0}. 

Another  observation  is  that  7r(Si)  is  a  compact  subset  of  D\.  Thus  by 
corollary  3.2,  S\  is  a  compact  subset  of  £ei. 

We  can  make  the  same  arguments  for  the  case  when  i  =  2  to  prove  that 

52  is  a  compact  subset  of  Se2  —  { (<72  ,£>2)  1^2  =  0}. 

Hence  by  letting  c  <  mm{ci,  C2},  it  is  true  that 

Sm  C  Si  x  S2  C  £ei  x  £e2  -  {z\Ai  =  0  or  A2  =  0}  (58) 

Thus,  S m  is  a  compact  subset  of  £ei  x  ^e2  —  {z\Ai  =  0  or  A2  =  0}  ■ 

We  can  now  apply  LaSalle’s  invariance  principle  to  show  that  the  trajec¬ 
tory  of  the  closed  loop  system,  starting  within  Sm,  converges  to  the  maximal 
invariant  subset  of  Sm  where  u(t)  =  0  is  satisfied  for  all  t. 

Proposition  5.2  With  V ,c  and  ut  given  as  in  (44)  ,(51)  and  (49),  the 
trajectory  starting  from  point  (qio,Pio,Q2o,P2o)  which  satisfies 


V{qio,Pio,q2o,P2o)  <  c 


(59) 
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will  converge  to  the  set  where 

k  = 

Ai  = 

(Mi -A h)  = 

are  satisfied  for  i  =  1,2.  Proof  In  order  to  calculate  the  invariant  set, let 
u\  =  0.  When  sin{E\ )  /  0,  we  get 

(ii  -  h\  +  Ci«w(— — J2  T -)/i)xgi 

+Zi x  (Hi  -  Hdi  +  pism(— — ^ -)Hi) 

Y  y  i 

+  ((Hi  -  Hdi  +  pism(  — - - ~)Ai)xpi)xqi 

=  0  (61) 

Take  inner  products  on  both  sides  with  q\  (t)  to  get 

(h  x  (Hi  -  Ad i  +  Cism(— - T -)Hi))  •  qi  =  0  (62) 

This  is  equivalent  to 

{h  x  qi(t))  ■  {Ax  -  Adl  +  Ci sin{— - ^2  T -)Hi)  =  0  (63) 

Let 

B  =  Ax  -  Adi  +  Ci sin{— — ^2  if. ^)Hi  (64) 

Equation  (63)  means  B  is  perpendicular  to  the  vector  lx  x  qx{t)-  We  can 
see  that  vector  B  should  stay  in  the  plane  spanned  by  lx  and  qx(t)- 

However,  from  the  assumption  that  ux{t)  =  0,  we  know  that  Ax  and  lx 
are  constant  vectors.  The  time  varying  vector  B  will  sweep  a  line  segment 
passing  the  fixed  point  (Hi  —  H^i).  The  direction  of  this  line  segment  is 
aligned  with  A\.  So  vector  B  will  be  the  intersection  of  the  {li,qi)  plane 
and  this  line  segment.  Because  q\  (t)  is  sweeping  the  orbital  plane,  the  (lx,qi) 
plane  is  identical  at  t  and  t  +  kTx  where  k  is  an  integer  and  Ti  is  the  period 
of  the  first  satellite.  Since  the  line  segment  is  not  changed,  the  intersection 
points  in  these  cases  must  be  identical.  Thus  we  must  have 


Idi 

Adi 

(j>  (60) 


B{t)  =B{t  +  kTx) 


(65) 


Coordinated  orbit  transfer  for  satellite  clusters 


15 


Figure  1:  The  relationship  between  h,q\,  A\  —  A ^  and  B 


Without  lost  of  generality,  suppose  at  time  t,  E\(t),  E2(t)  G  (— e,  ir  +  e). 
Then  equation  (65)  requires  that 

=  Ci{Ei(t  +  fcTi))sm(Tl(*~l~fcTl'>~^2(*+fcTl'l~^)  (66) 

Let  k  =  1  in  equation  (66),  because  (i(E\(t))  =  (i(Ei(t  +  T\)),  the  first 
observation  we  make  is  that  the  two  satellites  must  have  the  same  period.  In 
fact,  suppose  at  time  to  equation  (66)  is  satisfied.  Then  at  time  to  +  T\ ,  since 
Til  (to)  =  £i(to  +  T\),T  i(to)  =  Ti(to  +  Ti)  and  all  the  angles  (anomalies)  are 
in  the  range  of  [0,  27 r),  we  must  have  T2(to)  =  T2(to  +  T\).  But 

T2(to)  -  T2(to  +  Tx)  =  M(M2(t0  +  TO  -  M2(f0))  (67) 

V  ad 

Then  T2(to)  =  ^2(^0  +  Ti)  will  be  satisfied  only  if  M2(to  +  Ti)  =  M2(fo). 
Thus  we  shall  have  T\  =  k\  T2  where  k\  is  a  positive  integer.  Remember  we 
can  apply  the  same  argument  to  the  second  satellite  to  get  I2  =  k\T\ .  Thus 
we  must  have  T\  =  T^.  Hence  on  the  invariant  set,  we  proved  that  a\  =  0,2- 

On  the  other  hand,  for  a  specific  time  t,  we  know  that  there  exists 
t*  G  [0,  Ti)  such  that 


7T  +  /l(*)  =  + 


(68) 
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where  f±  is  the  true  anomaly  of  the  first  satellite.  The  value  of  t*  depends 
on  t.  The  plane  spanned  by  (h,qi)  at  time  f  will  also  be  identical  to  the 
plane  spanned  by  (h,qi)  at  time  t  +  t*.  Thus  we  must  have 


B(t)  =  B(t  +  t*) 

(69) 

which  requires  that 

Ci(JBi(t))Sm(Tl(t)-J2(t)-^) 

=  (i{Ei{t  +  t*))sin(ri{t+n~^2{t+n+4’) 

(70) 

Further,  ai  =  a 2  implies  that  M\(t)  —  M2(t)  =  M\{t  +  t *) 
one  can  verify  that 

-M2(t  +  t*), 

sin(r^t+r^t+n+<p) 

(71) 

For  (70)  to  be  satisfied,  one  possibility  is  that 

.  ,T!(t)-T2(*)-^  n 
sm{  2  )  =  0 

(72) 

Another  possibility  is  that 

Ci(E1(t))  =  -Ci(E1(t  +  t*)) 

(73) 

By  the  definition  of  (j,  one  can  verify  that  (73)  can  only  be  satisfied  when 
t  takes  value  from  a  set  of  measure  0.  Thus,  for  (70)  to  be  satisfied,  (72) 
must  be  true. 

Because  of  (72),  the  time  varying  parts  in  equation  (61)  vanish.  We  can 
make  the  same  argument  as  in  the  proof  of  the  single  satellite  case  [1]  to 
show  that 

h  =  Idl 

Tli  =  Adi  (74) 

We  can  apply  similar  arguments  for  the  second  satellite. 

Thus  we  have 

h  —  Idi 

- 1 ,  =  Adi 

(Tr-T  2)  =  ±0  (75) 
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for  i  =  1,2.  By  the  definition  of  Ti  and  T2  in  equation(38),  we  have 


(76) 


But  we  already  know  a\  =  (12  =  ad,  so  we  conclude  that 


(Mi  -M2)  =  <j> 


(77) 


6  Simulation  results 

To  verify  our  algorithm,  a  series  of  simulations  have  been  carried  out.  Here 
we  will  show  a  controlled  transfer  of  two  satellites  from  orbit  [a,  e,  i,  u ;,  0]  = 
[20, 0.1,  vr/4,  7t/2,  0]  with  initial  separation  of  mean  anomaly  being  7r/90  to 
the  orbit  [a,e,i,ui,Ci]  =  [25, 0.05,  7t/3,  vt/2,  0]  with  final  separation  of  mean 
anomaly  being  7t/18.  Only  relative  motion  between  the  satellites  are  plot¬ 
ted.  Figure  2  displays  the  desired  relative  motion  between  the  satellites. 
Figure  3  displays  the  relative  motion  between  the  satellites  using  the  con¬ 
trol  algorithm  proposed.  As  we  can  see,  the  desired  orbit  and  separation  are 
achieved. 


7  Future  work 

We  have  shown  that  the  shape  space  formed  from  the  angular  momentum 
vectors  and  Laplace  vectors  is  appropriate  to  describe  satellite  formations. 
The  control  laws  we  propose  are  based  on  a  Lyapunov  function  on  this 
shape  space.  It  has  several  limitations.  First,  we  have  not  consider  the 
case  when  the  formations  are  on  circular  orbits.  This  is  because  in  this  case 
the  Laplace  vectors  vanishes.  Second,  in  our  simulation  we  found  that  we 
can  drop  the  sin2(E)  terms  in  the  control  law  and  use  switching  laws  to 
improve  performance.  But  a  proof  of  the  convergence  requires  a  LaSalle’s 
invariance  principle  for  switching  laws  which  we  do  not  pursue  here.  Third, 
perturbations  such  as  J2  should  be  put  into  the  picture.  These  will  be 
investigated  in  future  work. 
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